library(SDMTools)
wd='/home/jc148322/rob/inputs/';setwd(wd)

#read in the asciis & create base.asc
bdall=read.asc('Bdall.asc')
bdwet=read.asc('Bdwet.asc')
bc12=read.asc('bc12.asc')
bc01=read.asc('bc01.asc')


#read in the data
pos=read.csv('coord.csv',as.is=TRUE)
threshold=read.csv('maxentResults.csv')
threshold=threshold$Minimum.training.presence.logistic.threshold
pos.frogs=read.csv('frog_pos.csv', as.is=TRUE)

########
##################################################################################################################
#CREATE THE LAYERS, EDIT THE DATA
out.dir='/home/jc148322/rob/outputs';setwd(out.dir)
#create the sig.diff ascii
bdall[which(bdall<threshold[1])]=0; #set values below threshold to 0
bdwet[which(bdwet<threshold[2])]=0; #set values below threshold to 0
base.asc=bdall; base.asc[which(base.asc>0)]=1

######################################
#calculate areas

bdall.binary=bdall; bdall.binary[which(bdall.binary>0)]=1;
cs.bdall=ClassStat(bdall.binary,cellsize=80, latlon=FALSE)
bdwet.binary=bdwet; bdwet.binary[which(bdwet.binary>0)]=1;
cs.bdwet=ClassStat(bdwet.binary,cellsize=80, latlon=FALSE)

area.names=c('original','new')
area.totals=c(cs.bdwet$total.area[2]/1000000,cs.bdall$total.area[2]/1000000)
areas=as.data.frame(area.totals,area.names)
write.csv(areas,'area.totals.csv')
######################################
#subset frog occurrence data and cbind bc01 and bc12
pos.frogs$bc01= extract.data(cbind(pos.frogs$Easting,pos.frogs$Northing),bc01)
pos.frogs$bc12 = extract.data(cbind(pos.frogs$Easting,pos.frogs$Northing),bc12)

frogdata=pos.frogs[,c('Species','bc01','bc12')]
lorica=frogdata[which(frogdata[,1]=='LITLORI'),];
nannotis=frogdata[which(frogdata[,1]=='LITNANN'),]; nannotis=nannotis[,2:3]

record.time=c('New','New','New','New','New','New','Old','Old','Old','Old','Old','Old','Old','Old','Old')
lorica=cbind(lorica,record.time)

######################################
#edit and simplify the data for biplot
tdata = cbind(pos$bc01,pos$bc12); #
#tdata[,1] = round(tdata[,1],1);
#tdata[,2] = round(tdata[,2]/10)*10 #round bc01 and bc12
tdata = unique(tdata) #find only unique combinations of rounded bc01 and bc12

pos$bdall[which(pos$bdall<threshold[1])]=0; #set values below threshold to 0
pos$bdwet[which(pos$bdwet<threshold[2])]=0;

wet.points=cbind(pos$bc01[which(pos$bdwet>0)],pos$bc12[which(pos$bdwet>0)])
all.points=cbind(pos$bc01[which(pos$bdall>0)],pos$bc12[which(pos$bdall>0)])
		
##################################################################################################################
#CREATE PLOTS AND IMAGES
plotcols=c('grey','blue','green')
#make the biplot
png('compare_frogs.png', width=8,height=8,units='cm', res=300,pointsize=5)

        plot(tdata, xlab="Mean Annual Temperature",ylab='Annual rainfall (mL)',xlim=range(tdata[,1],na.rm=T),ylim=range(tdata[,2],na.rm=T), type='n',cex.lab=1, cex.axis=1, font.lab='2')
        points(tdata, col='grey', pch=19,cex=3)

        points(all.points,col='blue', pch=19,cex=3)

        points(wet.points,col='green',pch=19,cex=3)

        points(nannotis,col='grey20',pch=4,cex=0.7,lwd=0.5)
				
        points(lorica$bc01[which(lorica$record.time=='New')],lorica$bc12[which(lorica$record.time=='New')],col='black',bg='orange',lwd=0.5, pch=21)
		points(lorica$bc01[which(lorica$record.time=='Old')],lorica$bc12[which(lorica$record.time=='Old')],col='black',bg='red',lwd=0.5,pch=21)
		
		legend(23,7000, c('Climate space', 'New Bd SDM', 'Original Bd SDM'), fill=plotcols, cex=1, bty='n')
		legend(23.1,6000, c('L. nannotis', 'L. lorica (old records)', 'L. lorica (new records)'), col=c('gray20','black','black'),pt.bg=c('orange','red'),pch=c(4,21,21), pt.lwd=0.5, cex=1, bty='n')
		
dev.off()

#nannotis points within wet and all
		poly.wet=chull(pos$bc01[which(pos$bdwet>0)],pos$bc12[which(pos$bdwet>0)])
		poly.wet=c(poly.wet,poly.wet[1])
		poly.all=chull(pos$bc01[which(pos$bdall>0)],pos$bc12[which(pos$bdall>0)])
		poly.all=c(poly.all,poly.all[1])
		
		nann.wet=pnt.in.poly(nannotis,wet.points[poly.wet,])
		nann.all=pnt.in.poly(nannotis,all.points[poly.all,])
		nann.dry=sum(nann.all$pip)-sum(nann.wet$pip)
nann.count=as.data.frame(rbind(c('nannotis.all','nannotis.wet','nannotis.dry'),c(sum(nann.all$pip),sum(nann.wet$pip),nann.dry)))
write.csv(nann.count,'nannotis.count.csv', row.names=FALSE)
######################################
#make the image

Colormap=c("grey",colorRampPalette(c('blue','deepskyblue','green','yellow','red'))(100))


pnts=cbind(x=c(470000,500000,500000,470000),y=c(8200000,8200000,8130000,8130000))

png ('maps2.png',width=2*dim(bdall)[1],height=1*dim(bdall)[2], bg= "white", pointsize=100)
    par(mfrow=c(1,2))
	
		image(bdwet, col=Colormap, ann=FALSE, axes=FALSE ,xlab="", ylab="")
		text(510000,8270000,'B',cex='1.5')
        image(bdall, col=Colormap,ann=FALSE, axes=FALSE ,xlab="", ylab="")
        legend.gradient(pnts,col=Colormap, c(0,1), title='Suitability',cex=1.3)
		text(510000,8270000,'C',cex='1.5')
dev.off()

#######################################
#plot frogpoints
nann.points=pos.frogs[which(pos.frogs$Species=='LITNANN'),]; nann.points=nann.points[,2:3]
lori.points=pos.frogs[which(pos.frogs$Species=='LITLORI'),]; lori.points=lori.points[,2:3]

png ('frog_points.png',width=1*dim(bdall)[1],height=1*dim(bdall)[2], bg= "white", pointsize=100)
    
        image(bdall, col=Colormap, ,xlab="", ylab="")
		points(nann.points[,2], nann.points[,1],col='grey20',pch=16, cex=0.7)
		points(lori.points[,2],lori.points[,1],  col='red',pch=16,cex=0.7)
		
        legend.gradient(pnts,col=Colormap, c(0,1), title='New',cex=1)
	
dev.off()



























